
x = [-1 1 1 -1];
y = [1 1 -1 -1];

vert = [x', y'];
edge = zeros(4, 2);
edge(:, 1) = [1:4]';
edge(:, 2) = mod(edge(:, 1), 4) + 1; 

param = createParam(0.00001, 0.03, 0.0000025, 0.1);

domain2D = createDomain(vert, edge, param);

dim = size(domain2D.vx, 1);

b = assemblyBoundaryVector(domain2D);

[rx, ry] = genInteriorNode(x, y, 10);

% res = 50;
% [rx, ry] = meshgrid(0:res,0:res);
% rx = rx / res * 1.99 - 0.995;
% ry = ry / res * 1.99 - 0.995;
% rx = reshape(rx, (res+1)*(res+1), 1);
% ry = reshape(ry, (res+1)*(res+1), 1);
% 
% [A, Z] = assemblyMatrix(domain2D);
% 
% theta = A \ b;
% 
% u = evaluateInteriorPoint([domain2D.vx; rx],[domain2D.vy; ry], domain2D, theta);
% 
% T = evaluateInteriorStress([domain2D.vx; rx],[domain2D.vy; ry], domain2D, theta);